week 7¶
In [1]:
import warnings
warnings.filterwarnings('ignore')
In [2]:
! jupyter nbextension enable --py --sys-prefix widgetsnbextension
In [3]:
import numpy as np
import pandas as pd
import datetime as dt
import scipy.stats as stats
In [4]:
# import holoviews as hv
# hv.extension('bokeh', logo=False)
# from holoviews.operation.datashader import aggregate
# import geoviews as gv
# import cartopy.crs as ccrs
In [5]:
from bokeh.plotting import figure, output_notebook, show
from bokeh.tile_providers import STAMEN_TERRAIN
output_notebook()
import datashader as ds
import datashader.transfer_functions as tf
from datashader.colors import Hot, viridis
from datashader.utils import export_image
from datashader.bokeh_ext import InteractiveImage
from holoviews.operation.datashader import aggregate
from functools import partial
from IPython.core.display import HTML, display
%matplotlib inline
import matplotlib.pyplot as plt
In [6]:
def process_file(file_name):
print('process of file: %s' %file_name)
# читаем исходный файл с данными
df = pd.read_csv(file_name, sep=',')
# в некоторые годы с 2014 столбцы назывались иначе и был пробел перед названием колонки, унифицируем
df.columns = [x.strip() for x in df.columns]
if 'vendor_id' in df.columns:
df.rename(columns={'vendor_id':'VendorID'}, inplace=True)
if 'pickup_datetime' in df.columns:
df.rename(columns={'pickup_datetime':'tpep_pickup_datetime'}, inplace=True)
if 'dropoff_datetime' in df.columns:
df.rename(columns={'dropoff_datetime':'tpep_dropoff_datetime'}, inplace=True)
if 'rate_code' in df.columns:
df.rename(columns={'rate_code':'RatecodeID'}, inplace=True)
if 'RateCodeID' in df.columns:
df.rename(columns={'RateCodeID':'RatecodeID'}, inplace=True)
# удаляем записи где длительность поездки равна 0
df = df[df.tpep_pickup_datetime != df.tpep_dropoff_datetime]
# удаляем записи где количество пассажиров равно 0
df = df[df.passenger_count != 0]
# удаляем записи где дальность поездки равна 0
df = df[df.trip_distance != 0.0]
# удаляем записи, не входящих в прямоугольник Нью-Йорка
df = df[(df.pickup_longitude <= -73.70001) & (df.pickup_longitude >= -74.25559) ]
df = df[(df.pickup_latitude >= 40.49612) & (df.pickup_latitude <= 40.91553) ]
# переводим время посадки пассажиров из формата object в datetime
df['tpep_pickup_datetime'] = pd.to_datetime(df['tpep_pickup_datetime'])
# переводим время посадки пассажиров из формата object в datetime
df['tpep_dropoff_datetime'] = pd.to_datetime(df['tpep_dropoff_datetime'])
# получаем время поездки для записи
df['trip_time'] = (df['tpep_dropoff_datetime'] - df['tpep_pickup_datetime']).dt.total_seconds()
# удаляем минуты или секунды из времени посадки пассажиров
df['tpep_pickup_datetime'] = df.tpep_pickup_datetime.values.astype('<M8[h]')
df.drop(['VendorID', 'passenger_count', 'trip_distance', 'RatecodeID', 'store_and_fwd_flag', 'fare_amount',
'extra', 'mta_tax', 'tip_amount', 'tolls_amount', 'improvement_surcharge', 'total_amount',
'trip_time'], axis=1, inplace=True)
# берем из файла с идентификаторми зон данные по 2500 зонам
regions_df = pd.read_csv('regions.csv', sep=';')
# добавляем в таблицу данные по идентификаторам зон посадки
binsx = regions_df.west[::50].values
binsy = regions_df.south[:50].values
res_pickup = stats.binned_statistic_2d(df.pickup_longitude, df.pickup_latitude, values=df.tpep_pickup_datetime,
statistic='count', bins=[binsx, binsy], expand_binnumbers=True)
df['region_pickup'] = (res_pickup.binnumber[0] - 1) * 50 + res_pickup.binnumber[1]
# добавляем в таблицу данные по идентификаторам зон высадки
# binsx = regions_df.west[::50].values
# binsy = regions_df.south[:50].values
# res_dropoff = stats.binned_statistic_2d(df.dropoff_longitude, df.dropoff_latitude, values=df.tpep_pickup_datetime,
# statistic='count', bins=[binsx, binsy], expand_binnumbers=True)
# df['region_dropoff'] = (res_dropoff.binnumber[0] - 1) * 50 + res_dropoff.binnumber[1]
# получаем дополнительный признак - зона города (из 265), файлы с описанием должны быть в папке ./taxi_zones
# !!!(этот участок обрабывается по 15 мин из-за spatial join
# zones_df = gpd.read_file('./taxi_zones/taxi_zones.shp').to_crs({'init': 'epsg:4326'})
# zones_df.drop(['OBJECTID', 'Shape_Leng', 'Shape_Area', 'zone', 'borough'], axis=1, inplace=True)
# shply_points = [Point(x, y) for x, y in zip(df['pickup_longitude'], df['pickup_latitude'])]
# res = gpd.sjoin(gpd.GeoDataFrame(crs={'init':'epsg:4326'}, geometry=shply_points), zones_df, how='left', op='within')
# df['zone'] = res['LocationID'].values
# специальный ряд для таблиц где нужно посчитать общее количество
df['num'] = np.ones(len(df))
df = df.set_index('tpep_pickup_datetime', drop=False)
return df
In [7]:
# %%time
# df_orig = pd.read_csv('yellow_tripdata_2016-06.csv')
In [8]:
# df_main_to_save = df_main[['pickup_longitude', 'pickup_latitude', 'region_pickup']]
In [9]:
# df_main_to_save.to_csv('yellow_tripdata_2016-06_edited.csv')
In [10]:
%%time
df_main = pd.read_csv('yellow_tripdata_2016-06_edited.csv', parse_dates=['tpep_pickup_datetime'])
df_main.set_index('tpep_pickup_datetime', inplace=True)
In [12]:
# загружаем предсказания количества поездок в июне 2016
predictions = pd.read_csv('predictions_2016-06.csv')
loc = predictions['id'].str.split('_')
regions = []; dates = []; hours = []; forwards = []
for x in loc:
regions.append(x[0])
dates.append(x[1])
hours.append(x[2])
forwards.append(int(x[3]))
df = pd.DataFrame()
df['regions'] = regions
df['datetime']= [date+' '+hour for date, hour in zip(dates, hours)]
df['predicts'] = forwards
df['datetime'] = [pd.datetime.strptime(row, '%Y-%m-%d %H') for row in df['datetime']]
df['y'] = predictions['y']
df.set_index('datetime', inplace=True)
In [13]:
def get_day_predicts(df, forward, region=None, date=None, hour=None):
if region is not None:
df = df[df['regions'] == str(region)]
df = df[df['predicts'] == int(forward)]
df.is_copy = False
df.index = df.index + dt.timedelta(hours=1)
if date is not None:
df = df[date]
if hour is not None:
df = df[df.index.hour == hour]['y'].values
return df
In [14]:
def get_bootstrap_samples(df_loc, predictions):
regions = df.regions.unique().astype(int)
result = {}
for n_samples, region in zip(predictions, regions):
n_samples = int(round(n_samples, 0))
data = df_loc[df_loc['region_pickup'] == region]
if n_samples != 0:
index = np.random.choice(data.shape[0], size=int(n_samples))
row = data.iloc[index, [0,1]].values
result[region] = row
else:
result[region] = [0, 0]
pred_coords = pd.DataFrame()
for x in result.values():
y = pd.DataFrame(x)
pred_coords = pd.concat([pred_coords, y])
# print(pred_coords.columns)
pred_coords.columns = ['pickup_longitude', 'pickup_latitude']
pred_coords.index = list(range(len(pred_coords)))
return pred_coords
In [15]:
def get_real_df(date, hour, region=None):
df_real = df_main[['pickup_longitude', 'pickup_latitude', 'region_pickup']]
if date is not None:
df_real = df_real[date]
if hour is not None:
df_real = df_real[df_real.index.hour == hour]
df_real = df_real[df_real['region_pickup'].isin(df.regions.unique().astype(int)) ]
if region is not None:
df_real = df_real[df_real['region_pickup'] == region]
df_real.drop('region_pickup', axis=1, inplace=True)
df_real.index = list(range(len(df_real)))
return df_real
In [16]:
def latlng_to_meters(df, lat_name, lng_name):
lat = df[lat_name]
lng = df[lng_name]
origin_shift = 2 * np.pi * 6378137 / 2.0
mx = lng * origin_shift / 180.0
my = np.log(np.tan((90 + lat) * np.pi / 360.0)) / (np.pi / 180.0)
my = my * origin_shift / 180.0
df.loc[:, 'pickup_x'] = mx
df.loc[:, 'pickup_y'] = my
In [17]:
import param, paramnb
In [18]:
df_loc = pd.read_csv('time_region_102_pickup_count.csv')
df_loc = df_loc.set_index('datetime')
df_loc = df_loc['2016-06-01':]
In [19]:
class GetData(param.Parameterized):
hour = param.Integer(default=12, bounds=(0,24))
date = param.ObjectSelector(default='2016-06-01',
objects=np.unique([x.strftime('%Y-%m-%d') for x in df_loc.index.astype('<M8[D]')]),
check_on_set=True)
1) карты с визуализацией реального и прогнозируемого спроса на такси в выбираемый пользователем момент времени¶
In [20]:
paramnb.Widgets(GetData, button=True, next_n=4)
In [21]:
date = GetData.date
hour = GetData.hour
# получаем координаты реальных проездок в выбранный день и час
df_real = get_real_df(date, hour, region=None)
# получаем предсказанное количество поездок в выбранный день и час (из ранее полученной модели)
if (date == '2016-06-30') & (hour > 18):
prediction_nums = get_day_predicts(df=df, forward=1, region=None, date=date, hour=hour)
print(prediction_nums.shape)
prediction_nums2 = get_day_predicts(df=df, forward=6, region=None, date=date, hour=hour)
print(prediction_nums2.shape)
prediction_nums = pd.concat([pd.DataFrame(prediction_nums), pd.DataFrame(prediction_nums2[-5:])], axis=0)
print(prediction_nums.shape)
# переиндексируем таблицу
hours = pd.date_range('2016-06-01 00:00:00', periods=720, freq='H')
prediction_nums.index = hours
else:
prediction_nums = get_day_predicts(df=df, forward=1, region=None, date=date, hour=hour)
# бутстрапируем коордианты из реальной выборки поячеечно для предсказанного количеста поездок в каждой ячейке
df_pred = get_bootstrap_samples(df_main, prediction_nums)
# добавляем в файл данные в метрах для получения карт ниже
latlng_to_meters(df_real, 'pickup_latitude', 'pickup_longitude')
latlng_to_meters(df_pred, 'pickup_latitude', 'pickup_longitude')
In [31]:
# NY = x_range, y_range = ((df_real['pickup_x'].min(), df_real['pickup_x'].max()),
# (df_real['pickup_y'].min(), df_real['pickup_y'].max()))
NYC = x_range, y_range = ((-8242000,-8210000), (4965000,4990000))
plot_width = 900
plot_height = 900
def base_plot(tools='pan, box_zoom, wheel_zoom, save, reset',
plot_width=plot_width, plot_height=plot_height, **plot_args):
p = figure(tools=tools, plot_width=plot_width, plot_height=plot_height,
x_range=x_range, y_range=y_range, outline_line_color=None,
min_border=0, min_border_left=0, min_border_right=0,
min_border_top=0, min_border_bottom=0, **plot_args)
p.axis.visible = False
p.xgrid.grid_line_color = None
p.ygrid.grid_line_color = None
return p
def create_image_pred(x_range, y_range, w=plot_width, h=plot_height):
cvs = ds.Canvas(plot_width=w, plot_height=h, x_range=x_range, y_range=y_range)
agg = cvs.points(df_pred, 'pickup_x', 'pickup_y')
img = tf.shade(agg.where(agg>np.percentile(agg,90)), cmap=viridis, how='eq_hist')
return tf.dynspread(img, threshold=0.5, max_px=4)
def create_image_real(x_range, y_range, w=plot_width, h=plot_height):
cvs = ds.Canvas(plot_width=w, plot_height=h, x_range=x_range, y_range=y_range)
agg = cvs.points(df_real, 'pickup_x', 'pickup_y')
img = tf.shade(agg.where(agg>np.percentile(agg,90)), cmap=viridis, how='eq_hist')
return tf.dynspread(img, threshold=0.5, max_px=4)
In [32]:
# карта с визуализацией реального спроса
export = partial(export_image, export_path="export")
p = base_plot()
p.add_tile(STAMEN_TERRAIN)
export(create_image_real(*NY), "NY1")
InteractiveImage(p, create_image_real)
Out[32]:
In [33]:
# карта с визуализацией предсказанного спроса
export = partial(export_image, export_path="export")
p = base_plot()
p.add_tile(STAMEN_TERRAIN)
export(create_image_pred(*NY), "NY2")
InteractiveImage(p, create_image_pred)
Out[33]:
2) временной ряд фактического и прогнозируемого спроса на такси в выбираемой области.¶
In [25]:
class GetCell(param.Parameterized):
cell = param.ObjectSelector(default=df.regions.unique().astype(int)[0],
objects=df.regions.unique().astype(int), check_on_set=True)
In [26]:
from sklearn.metrics import mean_absolute_error
In [27]:
paramnb.Widgets(GetCell, button=True, next_n=4)
In [34]:
region = GetCell.cell
# получаем предсказанное количество поездок в выбранный день и час (из ранее полученной модели)
prediction_nums = get_day_predicts(df=df, forward=1, region=region, date=None, hour=None)
# т.к. в файле предсказания на 6 часов вперед, то последние 5 записей не вошли в вышеукзанный файл, добавляем их
prediction_nums2 = get_day_predicts(df=df, forward=6, region=region, date=None, hour=None)
prediction_nums = pd.concat([prediction_nums, prediction_nums2[-5:]])
# переиндексируем таблицу
hours = pd.date_range('2016-06-01 00:00:00', periods=720, freq='H')
prediction_nums.index = hours
pred_arr = prediction_nums['y'].values
In [35]:
df_loc = pd.read_csv('time_region_102_pickup_count.csv')
df_loc = df_loc.set_index('datetime')
df_loc = df_loc['2016-06-01':]
real_arr = df_loc.loc[:,str(region)].values
In [36]:
plt.figure(figsize=(18,6))
plt.plot(pred_arr, 'r', label='predicted')
plt.plot(real_arr, 'b', label='real values')
plt.title('2016 June, cell: {}, MAE: {}, XGBRegressor'.format(region, mean_absolute_error(real_arr, pred_arr)))
plt.ylabel('pickups')
plt.xlabel('hours in june 2016')
plt.legend()
plt.show()
In [ ]: